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ABSTRACT 

Model-based prognostics approaches capture system knowl- 
edge in the form of physics-based models of components that 
include how they fail. These methods consist of a damage 
estimation phase, in which the health state of a component is 
estimated, and a prediction phase, in which the health state is 
projected forward in time to determine end of life. However, 
the damage estimation problem is often multi-dimensional 
and computationally intensive. We propose a model decom- 
position approach adapted from the diagnosis community, 
called possible conflicts, in order to both improve the com- 
putational efficiency of damage estimation, and formulate a 
damage estimation approach that is inherently distributed. 
Local state estimates are combined into a global state esti- 
mate from which prediction is performed. Using a centrifugal 
pump as a case study, we perform a number of simulation- 
based experiments to demonstrate the approach. 

1. Introduction 

Model-based prognostics approaches capture knowledge of 
how a system and its components fail through the use of 
physics-based models that capture the underlying physical 
phenomena (Daigle & Goebel, 2010b; Saha & Goebel, 2009; 
Luo, Pattipati, Qiao, & Chigusa, 2008). Model-based prog- 
nostics algorithms consist of two parts: (i) damage estima- 
tion, which is fundamentally a joint state-parameter estima- 
tion problem, and (ii) prediction, which projects the current 
joint state-parameter estimate forward in time to determine 
end of life (EOL). In (Daigle & Goebel, 201 1), we developed 
a prognostics framework using particle filters for the damage 
estimation step that handles several simultaneously progress- 
ing damage processes. However, the approach may not scale 
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well as the number of damage processes to track (i.e., the di- 
mension of the system) increases. 

In this paper, we improve both the scalability and the com- 
putational efficiency of the damage estimation task by ex- 
ploiting structural model decomposition (Williams & Mil- 
lar, 1998), similar to methods developed within the diagno- 
sis community (Bregon, Pulido, & Biswas, 2009; Roychoud- 
hury, Biswas, & Koutsoukos, 2009; Staroswiecki & Declerck, 
1989). In particular, we adopt the possible conflicts (PCs) 
approach (Pulido & Alonso-Gonzalez, 2004). PCs decom- 
pose a global system model into minimal overdetermined 
subsystems (local submodels) for fault detection and isola- 
tion (Pulido & Alonso-Gonzalez, 2004). PCs have also been 
used to formulate smaller estimation tasks for fault identifi- 
cation (Bregon, Pulido, & Biswas, 2009). In general, PCs 
can be used to automatically decompose a global joint state- 
parameter estimation task into a set of local estimation tasks 
that are easier to solve and require less overall computation. 
We use the PC approach to derive a minimal set of submod- 
els and define a local damage estimation task for each one. 
Every local estimator computes a local joint state-parameter 
estimate, represented as a probability distribution. Then, the 
local estimates are merged into a global estimate from which 
prediction is performed in the typical way (Daigle & Goebel, 
2011). 

The models are decomposed into independent submodels by 
using measured signals as local inputs. Therefore, each local 
estimator operates independently, and the damage estimation 
becomes naturally distributed. Clearly then, this approach es- 
tablishes a formal basis for distributed prognostics. This is in 
contrast to other proposed distributed prognostics approaches, 
e.g. (Saha, Saha, & Goebel, 2009), which still treat the prog- 
nostics problem as a global one in which only the computation 
is distributed, whereas we propose to decompose the global 
problem into a set of local ones for which computation may 
be trivially distributed. 
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Figure 1 . Prognostics architecture. 


We demonstrate our prognostics methodology on a centrifu- 
gal pump. Centrifugal pumps appear in a variety of domains 
and often need to be operated for long time periods, hence di- 
agnostics and prognostics become critical to ensuring contin- 
ued operation that meets performance requirements. We ap- 
ply our model-based prognostic approach based on structural 
model decomposition to centrifugal pumps using a number of 
simulation-based experiments when multiple damage mecha- 
nisms are active, and compare to results using the global esti- 
mation approach presented in (Daigle & Goebel, 201 1). 

The paper is organized as follows. Section 2 formally defines 
the prognostics problem and describes the prognostics archi- 
tecture. Section 3 describes the modeling methodology and 
develops the centrifugal pump model for prognostics. Sec- 
tion 4 presents the model decomposition approach and pro- 
vides results for the pump model. Section 5 describes the par- 
ticle filter-based local damage estimation method. Section 6 
discusses the prediction methodology. Section 7 provides re- 
sults from simulation-based experiments and evaluates the ap- 
proach. Section 8 concludes the paper. 

2. Prognostics Approach 

The goal of prognostics is the prediction of EOL and/or re- 
maining useful life (RUL) of a component. In this section, 
we first formally define the problem of prognostics. We then 
describe the model-based prognostics architecture based on 
structural model decomposition. 

2.1 Problem Formulation 

In general, we define a system model as 

x(f) = f(t,x(t),0(t),u(f),v(f)) 
y (t) = h(t, x(t), 0(t), u(t), n(t)), 

where t £ R is a continuous time variable, x(i) £ R™* is the 
state vector, 0(f) £ R" e is the parameter vector, u(f) € R" u 


is the input vector, v(f) £ is the process noise vec- 
tor, f is the state equation, y(f) £ R"» is the output vector, 
n(f) £ R n ™ is the measurement noise vector, and h is the 
output equation. The parameters 0(f) evolve in an unknown 
way. 

The goal is to predict EOL (and/or RUL) at a given time 
point tp using the discrete sequence of observations up to 
time tp, denoted as yo-.t P ■ The component must meet a 
given set of functional requirements. We say the compo- 
nent has failed when it no longer meets one of these require- 
ments. In general, we may capture this boundary on accept- 
able component behavior using a threshold that is a func- 
tion of the system state and parameters, TEOL( x (f), 0(f)), 
where TEoz,(x(f), 0(f)) = 1 if the system has failed and 
T E OL(x(t),0(t)) = 0 otherwise. Using T E ql , we formally 
define EOL with 

EOL{tp) = inf{f £ R : t > tp A TpoLip^it), 0(f)) = 1}, 

i.e., EOL is the earliest time point at which the damage thresh- 
old is met. RUL is then 

RUL{t P ) = EOL(t P ) - t P . 

Due to the noise terms v(f) and n(t), and uncertainty 
in the future inputs of the system, we at best compute 
only a probability distribution of the EOL or RUL, i.e., 

p(EOL(t P ) |yo: tp ) oTp(RUL(tp)\y 0 :t P ). 

2.2 Prognostics Architecture 

In our model-based approach, we develop detailed physics- 
based models of components and systems that include de- 
scriptions of how faults and damage evolves in time. These 
models depend on unknown parameters 0(t). Therefore, 
damage estimation is fundamentally a joint state-parameter 
estimation problem. In discrete time k, we jointly estimate 
x/, and Ok, and use these estimates to predict EOL and RUL 
at desired time points. Here, we assume that prognostics is 
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Figure 2. Centrifugal pump. 


not aided by a fault diagnosis module, and so we must jointly 
estimate all possible damage modes. 

We employ the prognostics architecture in Fig. 1. The sys- 
tem is provided with inputs u/, : and provides measured out- 
puts y/ c . The damage estimation module takes as input both 
u k and yk, and produces the estimate p(x k , 9k\yo:k )■ In this 
work, we decompose the global damage estimation problem 
into several subproblems based on model decomposition, as 
shown in Fig. 1. A model decomposition algorithm splits the 
global model into n submodels, Adi, A4 2 , ■ ■ ■ , A4 n . We con- 
struct for each submodel Ad,: a local estimator that performs 
damage estimation. Each estimator has input u). C rp, and 
y k ^ y k and produces the local state estimate p(xj.,0^.|yg. fc ), 
where xj. C x/, : and 0 ). C 6 k . Note that for two submod- 
els Adi and Adj, in general it is possible that xj. Pi xi ^ 0 
and 9 k Cl 6 k A 0. The local estimates are merged into the 
global estimate p(x kl 6 k ]yo:k)- The prediction module uses 
this joint state-parameter distribution, along with hypothe- 
sized future inputs, to compute EOL and RUL as probabil- 
ity distributions p(EOL kp \y 0:kp ) and p(RUL kp \y 0:kp ) at 
given prediction times kp. 

3. Pump Modeling 

We apply our prognostics approach to a centrifugal pump, and 
develop a physics-based model of its nominal and faulty be- 
havior. Centrifugal pumps are used in a variety of domains 
for fluid delivery. A schematic of a typical centrifugal pump 
is shown in Fig. 2. Fluid enters the inlet, and the rotation of 
the impeller, driven by an electric motor, forces fluid through 
the outlet. Radial and thrust bearings, along with lubricating 
oil contained within the bearing housing, helps to minimize 
friction along the pump shaft. Wear rings prevent internal 
pump leakage from the outlet to the inlet side of the impeller, 
but a small clearance is typically allowed to minimize fric- 
tion fa small internal leakage is normal). We review here the 
main features of the model, and refer the reader to (Daigle & 
Goebel, 2011) for details. 

The state of the pump is given by 

x(i)=[w(i) T t (t) T r (t) T 0 (t )] T , 


where u(t) is the rotational velocity of the pump, T t (t) is the 
thrust bearing temperature, T r (t ) is the radial bearing temper- 
ature, and T 0 (t) is the oil temperature. 

The rotational velocity of the pump is described using a 
torque balance, 

w = j (r e (f) - ru(t) - r L (t )) , 

where J is the lumped motor/pump inertia, r e is the electro- 
magnetic torque provided by the motor, r is the lumped fric- 
tion parameter, and r k is the load torque. 

We assume the pump is driven by an induction motor with 
a polyphase supply. A torque is produced on the rotor only 
when there is a difference, i.e., a slip, between the syn- 
chronous speed of the supply voltage, co H and the mechanical 
rotation, w. Slip, s, is defined as 

U) s — id 

S = . 

LO s 

The expression for the torque r e is derived from an equiva- 
lent circuit representation for the three-phase induction mo- 
tor, based on rotor and stator resistances and inductances and 
the slip s (Lyshevski, 1999): 

^ TipR 2 Irms 

sco s (f?i + R 2 /s) 2 + (uJsLx + u s L 2 ) 2 ’ 

where U \ is the stator resistance, L\ is the stator inductance, 
R 2 is the rotor resistance, L 2 is the rotor inductance, n is the 
number of phases, and p is the number of magnetic pole pairs. 
The dependence of torque on slip creates a feedback loop that 
causes the rotor to follow the rotation of the magnetic field. 
The rotor speed may be controlled by changing the input fre- 
quency oj s . 

The load torque 77 , is a polynomial function of the flow 
rate through the pump and the impeller rotational veloc- 
ity (Kallespe, 2005): 

tl = a 0 uj 2 + a k ujQ — a 2 Q 2 , 

where Q is the flow, and 00 , a\, and ci 2 are coefficients derived 
from the pump geometry (Kallespe, 2005). 

The rotation of the impeller creates a pressure difference from 
the inlet to the outlet of the pump, which drives the pump flow, 
Q. The pump pressure is computed as 

p p = Au> 2 + biuiQ - b 2 Q 2 , 

where A is the impeller area, and b\ and h> are coefficients 
derived from the pump geometry. Flow through the impeller, 
Qi , is computed using the pressure differences: 

Qi = C\j\p s +Pp- Pd\sign(p a + p p - Pd ), 

where c is a flow coefficient, p s is the suction pressure, and 
Pd is the discharge pressure. The small (normal) leakage flow 
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from the discharge end to the suction end due to the clearance 
between the wear rings and the impeller is described by 

Qi = ci\J\Vd -p s \sign(p d -p 3 ), 

where q is a flow coefficient. The discharge flow, Q, is then 

Q — Qi Qi ■ 

Pump temperatures are monitored as indicators of pump con- 
dition. The oil heats up due to the radial and thrust bearings 
and cools to the environment: 

To = ^-{HoATt - T 0 ) + H 0 ATr - T 0 ) 

Jo 

-H 0 i 3 (To-T a )), 

where J a is the thermal inertia of the oil, and the //„ , terms 
are heat transfer coefficients. The thrust bearings heat up due 
to the friction between the pump shaft and the bearings, and 
cool to the oil and the environment: 

T t = -'-(nut 2 - H t ATt - A) - H ta (T t - T a )), 

•>t 

where J t is the thermal inertia of the thrust bearings, r f is the 
friction coefficient for the thrust bearings, and the // t t terms 
are heat transfer coefficients. The radial bearings behave sim- 
ilarly: 

T r = — (ryw 2 — H r i(T r — T a ) — H r 2(T r — T a )), 

where the parameters here take on analogous definitions. 

The overall input vector u is given by 

u(f) = [p s (t) p d (t) T a (t) V(t) w s (f)] J . 

The measurement vector y is given by 

y(t) = [uj(t) Q(t) T t (t) T r (t ) T a (t)] T . 

Fig. 3 shows nominal pump operation. The input voltage (and 
frequency) are varied to control the pump speed. The electro- 
magnetic torque is produced initially as slip is 1. This causes 
a rotation of the motor to match the rotation of the magnetic 
field, with a small amount of slip remaining, depending on 
how large the load and friction torques are. As the pump ro- 
tates, fluid flow is created. The bearings heat up as the pump 
rotates and cool when the pump rotation slows. 

3.1 Damage Modeling 

For the purposes of prognostics, the model must include dam- 
age variables d C x representing the amount of particular 
forms of damage. The most significant forms of damage for 
pumps are impeller wear, caused by cavitation and erosion by 
the flow, and bearing failure, caused by friction-induced wear 
of the bearings. In each case, we map the damage to a par- 
ticular parameter in the nominal model, and this parameter 
becomes a state variable in d(f). The evolution of these dam- 
age variables is described by damage progression equations. 


Input Voltage 


Rotational Velocity 







Figure 3. Nominal pump operation. 

which are included in the state equation f. These equations 
are parameterized by unknown wear parameters w CO. 

Impeller wear is represented as a decrease in effective im- 
peller area A (Biswas & Mahadevan, 2007; Daigle & Goebel, 
201 1). We use the erosive wear equation (Hutchings, 1992): 

A{t) = - w a Q 2 , 

where w a is a wear coefficient. A decrease in the impeller 
area will decrease the pump pressure, which, in turn, reduces 
the delivered flow, and, therefore, pump efficiency. The pump 
must operate at a certain minimal efficiency, defining an EOL 
criteria. We define A~ as the minimum value of the impeller 
area at which this requirement is met, hence, Teol = 1 if 
A(t) < A- . 

Bearing wear is captured as an increase in friction. Sliding 
and rolling friction generate wear of material which increases 
the effective coefficient of friction (Hutchings, 1992; Daigle 
& Goebel, 2010b, 2011): 

h(t) = w t r t u 2 , 
f r (t) = w r r r uj 2 , 

where wt and w r are the wear coefficients. The slip com- 
pensation provided by the electromagnetic torque generation 
will mask small changes in friction, but these changes can be 
observed using the bearing temperatures. Limits on the max- 
imum values of these temperatures define EOL for bearing 
wear. We define r) f and r + as the maximum permissible val- 
ues of the friction coefficients, before the temperature limits 
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Figure 4. System graph and minimal subgraphs of the pump. 


are exceeded over a typical usage cycle. So, Teol = 1 if 
r t (f) > or r r (t) > r+. 

So, the damage variables are given by 

d(t) = [A(t) r t (t) r r {t)\ T , 

and the full state vector becomes 

x(f) = [w(t) T t (t) T r (t) T 0 (t) A(t) r t (t ) r r (t)] T . 

The wear parameters form the unknown parameter vector, i.e., 

w(f) = 0(t) = [wa w t w r ] T . 

4 . Model Decomposition 

Especially for nonlinear systems, state and parameter estima- 
tion is a nontrivial problem for which no general closed-form 
solution exists. In general, these problems may be solved 
by numerical optimization methods, where the computational 
complexity is exponential in the size of the model, or by 
approximate Bayesian filtering methods, like particle filters, 
where the computational complexity is only linear in the size 
of the sample space, but, the number of sufficient samples 
grows with the model size. 

Several approaches have been developed to decrease the com- 
putational complexity by model decomposition, in which the 
global model is decomposed into several independent sub- 
models (Staroswiecki & Declerck, 1989; Williams & Millar, 


1998). This results in a set of smaller, lower-dimensional es- 
timation tasks. We adopt the PC approach for model decom- 
position. PCs are minimal subsets of equations with sufficient 
analytical redundancy to generate fault hypotheses from ob- 
served measurement deviations, and, in previous work, we 
used PCs to propose a more robust and computationally sim- 
pler parameter estimation approach for fault identification 
(Bregon, Pulido, & Biswas, 2009), where the parameter es- 
timation task using the entire system model was replaced 
by a set of smaller estimation problems (one for each PC). 
In this approach, fault identification is fundamentally a joint 
state -parameter estimation problem. This is equivalent to the 
damage estimation problem in prognostics, only the models 
specifically include damage progression. So, in this paper, 
we adopt this paradigm for distributed damage estimation. 

In order to compute the minimal set of submodels of a sys- 
tem, a structural representation of the system model is needed. 
In previous work, we computed submodels as PCs using hy- 
pergraphs (Pulido & Alonso-Gonzalez, 2004) or Temporal 
Causal Graphs (TCGs) (Bregon, Pulido, & Biswas, 2009) as 
inputs. Representations that include computational causal- 
ity, like TCGs, are favored because causality allows efficient 
derivation of PCs. In this work, we represent the system 
model with a directed hypergraph, and use an algorithm close 
to the TCG-based algorithm presented in (Bregon, Pulido, 
Biswas, & Koutsoukos, 2009). 
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The system model M. is represented using a set of functions 
T over variables V, where a subset of the variables X C V 
corresponds to the state variables x, a subset 0 C V corre- 
sponds to the unknown parameters 6 , a subset U C V corre- 
sponds to the (known) inputs u, and a subset Y C V corre- 
sponds to the (measured) outputs y. For the purposes of the 
model decomposition algorithm, we represent M using an 
extended directed hypergraph Q = (V. E, F), where V is the 
set of vertices corresponding directly to the variables in Al, 
E is the set of hyperedges of the form (V 1 , v) with V' C V 
being a set of vertices and v £ V being a single vertex, and F 
is a map from an edge (V 7 , v) £ E to a function f £ T such 
that v = f(v !,V 2 , ■ ■ ■ ,v n ) where V' = {i’i,v 2 , ■ ■ -,v n }. M 
and Q are equivalent data structures, where the direction of the 
edges in Q captures the computational causality of the func- 
tions in T . The variable sets X , 0, Y, and U are represented 
equivalently in Q as vertex sets, and we use y, £ Y to refer to 
the vertex/variable corresponding to the ?'th output in y. 

Fig. 4a shows the hypergraph Q for the pump model described 
in Section 3. Individual arrows pointing to the same vertex are 
to be interpreted as a single directed hyperedge. State vari- 
ables are denoted using dashed circles, and measured vari- 
ables are denoted with boxes. 

Algorithm 1 computes subgraphs corresponding to PCs from 
a system graph Q. One PC will be computed for each output 
y t £ Y. So, for the pump model, there will be five total sub- 
models, one each for at, Q, T t , T r , and T 0 . For each vertex 
in Y, the algorithm propagates back to members of U and Y, 
which will be used as inputs to the submodel. Vertices which 
are not included in U or Y or have not yet been included in 
Vi are added to the set vertices for further backward propa- 
gation. For example, starting with T t in Fig. 4a to form the 
subgraph shown in Fig. 4d, we propagate to T a , a measured 
output, at which propagation terminates, and T t . From T t we 
propagate back further to the input T a and the state r t , from 
which we propagate to r t , from which we propagate to mea- 
sured output oj and the unknown parameter wt . All vertices 
and edges encountered are added to Vi and E t , respectively, 
and the model equations corresponding to the added edges are 
included in the submodel as well. The algorithm forms from 
Q subgraphs Gi — (Vi, Ei, Fi), which may be easily trans- 
lated to submodels AT;. If each v £ X U 0 i s causally linked 
to at least one output, then every variable x £ X and 9 £ 0 
will belong to at least one V t over the set of G, computed, i.e., 
will belong to at least one submodel AT;. 

The algorithm decomposes the pump model into 5 submodels, 
with their corresponding subgraphs shown as Figs. 4b to 4f. 
For example, the T t subgraph takes as input measurements of 
u> and T 0 , and computes the expected value of T t . Damage 
estimation for this submodel will compute estimates of T t , 
rt, and wt.. Note that, for the pump model, each state or pa- 
rameter will be estimated by exactly one submodel, therefore, 
there will be no overlap in the local estimates. 


Algorithm 1 {Gi}™Zi = Decompose^) 

for ( = 1 to n y do 

Vi {Vi} 

Ei £- 0 

vertices <— {yi} 
while vertices 0 do 
v vertices{ 1} 
vertices <— vertices \ {v} 
edges <- {( V',v ) £ E : V' C V} 
for all (V', v) £ edges do 
for all v' £ V' do 

if v' U and v' Y and v' <jL V then 

vertices <— vertices U {t/} 

end if 
end for 

Vi <- Vi U V' 

Ei<-EiU{(V',v)} 

Fi(V',v) ■£- F(V',v) 

end for 
end while 

Qi <- ( Vi,Ei,Fi ) 

end for 


5. Damage Estimation 

In our local estimation scheme, the local estimator for each 
submodel AT; produces a local estimate p(x^, |yo:fe), 
where xj, C x k and 0 l k C 6 k . The local estimates are com- 
bined into the global state estimate p(x k , 0 k |yo:fc)- 
Due to the decoupling introduced by the decomposition 
scheme, we lose information about the covariance between 
states in separate submodels. If these covariances are nom- 
inally small, then this information loss is acceptable and the 
approximation of the global state estimate obtained by merg- 
ing the local state estimates into a global state estimate will 
closely approximate the global estimate obtained through a 
global damage estimator. Although we lose information due 
to the decoupling, the advantage is that the local estimation 
tasks are naturally distributed, and therefore, unlike the global 
estimation approach, the distributed approach scales well as 
the size of the model increases. Further, the local estimation 
tasks become easier to solve and should require less compu- 
tational resources without sacrificing estimation performance. 
This should be the case as long as the sensor measurements 
that are used as inputs to the submodels are reliable and do 
not exhibit extremely high noise. 

A general solution to the problem of damage estimation is the 
particle filter, which may be directly applied to nonlinear sys- 
tems with non-Gaussian noise terms (Arulampalam, Maskell, 
Gordon, & Clapp, 2002). The main disadvantage of the par- 
ticle filter is the computational complexity, which is linear 
in the amount of samples, or particles, that are used to ap- 
proximate the state distribution, as typically a large number 
of particles are needed, and the sufficient number of particles 
increases with the dimension of the state-parameter space. A 
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key advantage of the model decomposition algorithm is that 
it creates submodels that are simpler than the global model. 
Some of these submodels may be completely linear, and some 
may require only state estimation, not joint state-parameter 
estimation. If only state estimation is required, the Kalman 
filter or one of its nonlinear extensions may be used, which 
requires computation on the order of a particle filter using 
one particle (e.g., Kalman filter or extended Kalman filter) or 
a number of particles linear in the state dimension (e.g., un- 
scented Kalman filter (Julier & Uhlmann, 1997 )), resulting in 
a significant improvement in computational complexity. The 
remaining submodels that require joint state-parameter esti- 
mation represent several small, low-dimensional estimation 
problems that are easier to solve than the global one, there- 
fore requiring less computation overall. 

In particle filters, the state distribution is approximated by a 
set of discrete weighted samples, or particles: 

where N denotes the number of particles, and for particle i, 

(i) (i) 

x k denotes the state vector estimate, 0 k denotes the pa- 

rameter vector estimate, and w k denotes the weight. The 
posterior density is approximated by 

N 

p(x k , 6 k |y 0: fc) « 0(Ox(dx fc d0fc), 

i— 1 

where 5, <->) o a),(dx k d0 k ) denotes the Dirac delta function 
located at (x£ 0 k ^). 

We use the sampling importance resampling (SIR) particle fil- 
ter, using systematic resampling. Each particle is propagated 
forward to time k by first sampling new parameter values, 
and then sampling new states using the model. The particle 
weight is assigned using y k . The weights are then normal- 
ized, followed by the resampling step. Pseudocode is pro- 
vided in (Arulampalam et al., 2002 ; Daigle & Goebel, 201 1 ). 
The parameters 0 k evolve by some unknown random process 
that is independent of the state x fc . To perform parameter esti- 
mation within a particle filter framework, we assign a random 
walk evolution, i.e., 8 k = 0 k -i+£ k -i, where £, k _ 1 is a noise 
vector. During the sampling step, particles are generated with 
parameter values that will be different from the current values 
of the parameters. The particles with parameter values closest 
to the true values should match the outputs better, and there- 
fore be assigned higher weight. Resampling will cause more 
particles to be generated with similar values, so the particle 
filter converges to the true values as the process is repeated 
over each step of the algorithm. In general, though, conver- 
gence is not always guaranteed. 

The selected variance of the random walk noise determines 
both the rate of this convergence and the estimation perfor- 
mance after convergence. Therefore, this parameter should 


be tuned to obtain the best possible performance, but the op- 
timal value is dependent on the value of the hidden wear 
parameter, which is unknown. We use the variance control 
method presented in (Daigle & Goebel, 2011). In this ap- 
proach, the variance of the hidden wear parameter estimate is 
controlled to a user-specified range by modifying the random 
walk noise variance. We assume that the £ values are tuned 
initially based on the maximum expected wear rates. The al- 
gorithm uses relative median absolute deviation (RMAD) as 
the measure of spread. The adaptation scheme controls the 
error between the actual RMAD of a parameter 8(j), denoted 
as Vj , and the desired RMAD value (e.g., 10%), denoted as 
v*, using a proportional control strategy governend by a gain 
P (e.g., 1 x 10 3 ). There are two different setpoints. The 
first allows for a convergence period, with setpoint v* 0 (e.g., 
50%). Once Vj reaches T (e.g., l.2v* 0 ), a new setpoint v*^ 
(e.g., 10%) is established. The advantage of this methodol- 
ogy is that the random walk variance is automatically tuned to 
achieve the best performance for the requested relative spread 
for the actual value of the hidden parameter. 

6. Prediction 

Prediction is initiated at a given time kp. In order to obtain a 
prediction that is valid for the global state -parameter vector, 
we must first combine the local estimates into a global esti- 
mate. To do this, we assume that the local and global state 
estimates may be sufficiently approximated by a multivari- 
ate normal distribution A /”(/r, S). So, for each local state- 
parameter distribution i, we obtain the mean //' and covari- 
ance matrix S\ We then combine all of these into a global 
mean fi and covariance S. If there is overlap in the state- 
parameter estimates, i.e., if two submodels both estimate the 
same state variable x or parameter 9, then we take the average 
value for common means and covariances (alternate strategies 
may also be used). The covariance information lost due to the 
decoupling will appear as zeros in the global covariance ma- 
trix. 

We then sample from the global state-parameter distribution 
defined by Af(n, £) a number of times to obtain a set of sam- 
ples that sufficiently approximates the distribution defined by 
those parameters, which may each be simulated to EOL. An 
alternate approach is to use the unscented transform to de- 
terministically select the minimal number of samples from 
this distribution that capture the statistical moments as de- 
scribed in (Daigle & Goebel, 2010a). Although the latter 
method is computationally more efficient, we use the former 
method here in order to obtain a fair comparison to the re- 
sults presented for the global estimation approach in (Daigle 
& Goebel, 2011). 

Using the global joint state-parameter estimate at kp, 
p(x kp ,8 kp \yo :kp ), which represents the current knowl- 
edge of the system at time kp, the goal is to compute 
p{EOL k p\y 0:k p) and p(RUL kp |y 0;fcp ). As discussed in 
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Algorithm 2 EOL Prediction 


Inputs: {(xj: 0 


/)(*) 
^kp 5 k p 

Outputs: {EOL^ p 
for i = 1 to N do 

k ^ — k p 




,(0 


,(0 
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EOL' 
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(i) 


Section 5, the particle filter computes 

N 

p(x kp ,Ok P \yo:kp) Z w kl S ( x (i) e (» ddx kp dG kp ). 

‘ ^ ' kp' kp' 

2=1 

We can approximate a prediction distribution n steps forward 
as (Doucet, Godsill, & Andrieu, 2000) 

Pi^-kp+m @ kp+n\yO:kp') ~ 

N 

^ ', w k c d /■ (») g(<) \{d~K. kp j rn dd kp j rn ). 

v fcp+n’ fcp+n/ 

2=1 

So, for a particle i propagated n steps forward without new 

. u\ 

data, we may take its weight as w kp . Similarly, we can ap- 
proximate the EOL as 

N 

p(EOL kp |y 0;fcp ) ~J2 w kl S EOL^ i.dEOL kp ). 
i = i kp 

To compute EOL, then, we simulate each particle forward 
to its own EOL and use that particle’s weight at kp for 
the weight of its EOL prediction. The pseudocode for the 
prediction procedure is given as Algorithm 2 (Daigle & 
Goebel, 2010b). Each particle i is propagated forward until 
Teol(*-u\ Q\^) evaluates to 1. The algorithm hypothesizes 
future inputs of the system, u;,. . In this work, we consider the 
situation where a single future input trajectory is known. 

7. Results 

We performed a number of simulation-based experiments to 
analyze the performance of the prognostics approach using 
local damage estimation. For the purposes of comparison, we 
include results from the global estimation approach. For the 
local estimation approach, we use the five submodels derived 


in Section 4, where the three submodels associated with dam- 
age models use particle filters for joint state-parameter esti- 
mation, and the remaining two submodels, which require only 
state estimation, use extended Kalman filters. The global ap- 
proach used N = 500, so when the three local particle filters 
each use N = 167, the total computational cost is equivalent 
to that of the global particle filter. We try also N = 100 and 
N = 50 to observe the changes in performance when less to- 
tal computation is performed. Since the local estimators use 
measured values as inputs, performance will degrade as sen- 
sor noise is increased. We varied the sensor noise variance by 
factors of 1, 10, 100, and 1000, to explore this situation. 

In a single experiment, combinations of wear parameter val- 
ues were selected randomly within a range. We selected the 
true wear parameter values in [1 x 10” 3 ,4 x 10 -3 ] for ui a, 
and in [1 x 10 -11 , 7 x 10~ n ] for wt and w r , such that the 
maximum wear rates corresponded to a minimum EOL of 20 
hours. The local estimators had to estimate both the local 
states and the local unknown wear parameters. In all exper- 
iments, we used T = 60%, Vq = 50%, = 10%, and 

P = 1 x 10 4 for the variance control algorithm. We per- 
formed 20 experiments for each value of N and sensor noise 
level. We considered the case where the future input of the 
pump is known, and it is always operated at a constant RPM, 
in order to limit the uncertainty to only that involved in the 
noise terms and that introduced by the filtering algorithms. 
The averaged estimation and prediction performance results 
are shown in Table 1. The part of the table with \M\ = 1 
corresponds to results using the global model. The column 
labeled N lists the number of particles used per submodel, 
and the column labeled n lists the sensor noise variance 
multipliers. Here, we use percent root mean square error 
(PRMSE) as a measure of estimation accuracy, relative ac- 
curacy (RA) (Saxena, Celaya, Saha, Saha, & Goebel, 2010) 
as a measure of prediction accuracy, and RMAD as a mea- 
sure of spread. Each are averaged over multiple prediction 
points for a single scenario (see (Saxena et ah, 2010; Daigle & 
Goebel, 2011) for the mathematical definitions of the metrics 
used here). Note that all metrics are expressed as percentages. 
We can see that in the case where N = 167, for the same 
amount of computation, the local estimation approach obtains 
results very close to the global approach for damage estima- 
tion accuracy. In some cases, performance is slightly bet- 
ter, and in other cases, slightly worse. At the highest noise 
level, the local approach improves significantly for estima- 
tion of tv a. This is mostly due to the convergence proper- 
ties of the global approach. It tends to converge with much 
more difficulty than the local approach for the same amount 
of noise. RMADs of the wear parameters are also quite evenly 
matched. Here again, in some cases the local approach is 
slightly better, and in others slightly worse. As the sensor 
noise increases, the variance is naturally larger and more dif- 
ficult to control, resulting in the increase in RMAD as sen- 
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Table 1 . Estimation and Prediction Performance 


\M\ 

N 

n 

PRMSEu,^ 

PRMSE Wt 

PRMSE„ Jr 

RMAD wa 

RMAD Wt 

RMAD Wr 

RA 

R M A D at: /\ 

1 

500 

i 

3.70 

3.58 

2.54 

11.58 

11.27 

10.03 

97.28 

11.61 

1 

500 

10 

4.15 

2.81 

2.74 

12.25 

11.48 

10.63 

96.58 

12.34 

1 

500 

100 

6.30 

3.46 

3.23 

13.46 

12.38 

11.59 

94.69 

14.09 

1 

500 

1000 

12.93 

6.25 

5.29 

13.92 

12.99 

12.64 

79.37 

15.32 

5 

167 

1 

3.19 

2.61 

2.88 

12.26 

10.85 

10.76 

96.61 

12.01 

5 

167 

10 

3.66 

2.90 

3.56 

12.69 

11.09 

11.85 

95.28 

13.32 

5 

167 

100 

4.44 

3.39 

3.78 

13.05 

11.78 

12.56 

93.17 

14.57 

5 

167 

1000 

5.59 

4.46 

8.26 

14.72 

12.86 

15.09 

87.66 

16.19 

5 

100 

1 

4.02 

3.49 

3.49 

12.42 

10.68 

10.77 

95.90 

12.15 

5 

100 

10 

4.52 

3.78 

4.27 

12.69 

11.04 

11.45 

95.16 

13.30 

5 

100 

100 

5.71 

4.02 

6.05 

12.99 

11.81 

12.73 

91.82 

14.74 

5 

100 

1000 

9.06 

4.76 

6.91 

13.83 

12.15 

13.92 

79.90 

16.40 

5 

50 

1 

5.66 

4.98 

5.19 

12.33 

10.41 

10.39 

94.59 

12.39 

5 

50 

10 

6.12 

4.92 

6.29 

12.41 

10.71 

11.21 

93.44 

12.99 

5 

50 

100 

7.43 

6.08 

8.24 

12.94 

11.16 

11.91 

90.05 

14.19 

5 

50 

1000 

14.03 

9.33 

14.41 

12.90 

11.66 

12.46 

73.05 

13.62 


sor noise increases for both approaches. As the number of 
particles used in the local approach is reduced, accuracy de- 
creases, as expected, but not significantly. The RMADs actu- 
ally generally decrease as the number of particles is reduced, 
corresponding to increased precision at the cost of decreased 
accuracy. 

Prediction performance corresponds to the change in dam- 
age estimation performance. More accurate damage estimates 
correspond to higher RA, and increases in spread of the wear 
parameters leads to increases of the spread of the RUL. For 
N = 167, the performance is slightly worse, but still com- 
parable to the global estimation approach, even though some 
of the covariance information in the global state estimate was 
lost due to the decoupling. At the highest noise level, the local 
approach has significantly better accuracy. This is due to the 
relatively poor convergence behavior of the global approach, 
leading to inaccurate early predictions that bring down the av- 
eraged RA. As N decreases, so that less total computation is 
being performed, accuracy reduces, but so does spread. For 
N = 100, less computation is performed with only a small 
change in performance. As N is reduced further to 50, per- 
formance begins to degrade. Moreover, as sensor noise in- 
creases, the local approach can lose its advantage over the 
global approach, and this occurs when N is reduced to 50. 
For N = 50 and the nominal amount of sensor noise, com- 
parable prognostics performance is achieved to the global ap- 
proach with less than a third of the computation. We expect 
the benefits of the local approach to be more pronounced as 


the dimension of the state-parameter space increases. 

8. Conclusions 

In this paper, we developed a novel distributed damage esti- 
mation approach for model-based prognostics that is based on 
a formal framework for structural model decomposition. Us- 
ing the concept of PCs, a system model is decomposed into a 
set of minimal submodels. A local damage estimation prob- 
lem is defined for each submodel. Focal state-parameter es- 
timates obtained using an appropriate filter are merged into 
a global state-parameter estimate from which EOF predic- 
tions are computed. Results demonstrate that equivalent, or in 
some cases, better prognostics performance can be achieved 
using this methodology with less computation than a global 
approach. Further, the approach can be naturally distributed 
and therefore may serve as a fundamental aspect of a practical 
system-level prognostics approach. 

The idea of using model decomposition to improve state 
and parameter estimation is not new. For example, sub- 
space methods (Katayama, 2005) for system identification 
employ QR-factorization and singular-value decomposition 
(Overschee & Moor, 1996) for solving identification prob- 
lems in large-dimension systems. These methods are nu- 
merically robust for linear systems. Recently, several ex- 
tensions have been proposed that apply to nonlinear systems 
(e.g., (Westwick & Verhaegen, 1996)). However, methods 
to automatically derive the decomposition from the system 
model have not been addressed. 
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An approach for decomposing a system model into smaller 
hierarchically organized subsystems, called dissents, is de- 
scribed in (Williams & Millar, 1998). PCs are conceptu- 
ally equivalent to dissents, and previous work applied PCs 
for model decomposition to generate a more robust and com- 
putationally simpler parameter estimation approach for fault 
identification (Bregon, Pulido, & Biswas, 2009). Simulation 
results in that case showed an improvement in estimation ac- 
curacy while having a faster convergence to true solutions. 
Similar work was proposed in (Roychoudhury et al., 2009) 
using a dynamic Bayesian network (DBN) modeling frame- 
work, in which an automatic approach for model decompo- 
sition into factors based on structural observability was de- 
veloped for efficient state estimation and fault identification. 
This approach also obtained an improvement in state esti- 
mation efficiency without compromising estimation accuracy. 
The relation between both approaches has been established in 
(Alonso-Gonzalez, Moya, & Biswas, 2010), where DBNs are 
derived from PCs for the purposes of estimation. 

In future work, we will investigate extensions to system-level 
and distributed prognostics. For one, the model decomposi- 
tion algorithm suggests a sensor placement strategy to opti- 
mize the decomposition of a system-level model into inde- 
pendent component models. One need only place sensors at 
the inputs and outputs of components to ensure that compo- 
nent models may be decoupled. The submodels derived from 
the PC approach are minimal in that, for a given output, they 
contain only the subset of the model required to compute that 
output as a function of only inputs and other measured out- 
puts. Nonminimal submodels may be formed by merging 
minimal submodels, and this may be desired in some cases, 
e.g., if it eliminates using some high-noise sensors as inputs. 
This forms part of a more generalized model decomposition 
framework under development. As described in the paper, 
distributed damage estimation is an essential part of a dis- 
tributed model-based prognostics architecture. The computa- 
tion associated with the prediction problem in our approach 
can be trivially distributed (via parallel EOL simulations), but 
in future work, we would like to develop a decomposition of 
the prediction problem into local prediction problems for a 
fully distributed prognostics architecture. 
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